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A hybrid method has been developed for simulations of compressible turbulent mixing 
layers. Such mixing layers dominate the flows in exhaust systems of modern day aircraft 
and also those of hypersonic vehicles currently under development. The hybrid method 
uses a Reynolds-averaged Navier-Stokes (RANS) procedure to calculate wall bounded 
regions entering a mixing section, and a Large Eddy Simulation (LES) procedure to 
calculate the mixing dominated regions. A numerical technique was developed to enable 
the use of the hybrid RANS-LES method on stretched, non-Cartesian grids. The hybrid 
RANS-LES method is applied to a benchmark compressible mixing layer experiment. 
Preliminary two-dimensional calculations are used to investigate the effects of axial grid 
density and boundary conditions. Actual LES calculations, performed in three spatial 
directions, indicated an initial vortex shedding followed by rapid transition to turbulence, 
which is in agreement with experimental observations. 


Introduction 

The use of computational fluid dynamics (CFD) to 
assist in the analysis and design of aerospace vehicles 
and their components has substantially increased in re- 
cent years. For analyzing one particular class of flows, 
that of aircraft engine exhaust nozzles, Reynolds- 
averaged Navier-Stokes (RANS) codes have been used 
extensively by government organizations (i.e. NASA) 
and aerospace companies. Exhaust nozzles being de- 
veloped for modern day commercial aircraft typically 
have multiple streams with a core flow and one or more 
bypass streams which mix with the high energy core 
flow before exiting the nozzle to lower jet noise while 
maintaining high thrust levels. Similarly in NASA’s 
recent High-Speed Research program, the engine ex- 
haust systems for this proposed supersonic transport 
were designed to be mixer-ejector nozzles, which en- 
train secondary air into the exhaust nozzle to mix with 
the core engine stream, again with the goal of simul- 
taneously lowering jet noise and maintaining sufficient 
thrust. 1 Propulsion systems currently under develop- 
ment for use on hypersonic and reusable space launch 
vehicles, such as the Turbine-Based Combined-Cycle 
(TBCC) and Rocket-Based Combined-Cycle (RBCC) 
concepts also employ mixer ducts. 

The flows in these nozzle systems all have compress- 
ible turbulent mixing as the dominant flow character- 
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istic. RANS codes used by research and development 
engineers to analyze these nozzles have employed tur- 
bulence models to replace the unsteady turbulent mo- 
tion with an effective eddy viscosity. Unfortunately, no 
turbulence model has been developed to date which is 
able to accurately represent the turbulent motion for 
such nozzle flows. Validation studies 2,3 have show T n 
that the “state of the art” turbulence models avail- 
able in production-use RANS codes have major defi- 
ciencies in predicting turbulent mixing in nozzle and 
jet flows involving compressibility, high temperatures, 
and three-dimensionality. 

The known limitations of RANS techniques to cal- 
culate complex turbulent flows, coupled with contin- 
ually increasing computing power, have led to inter- 
est in more sophisticated calculation techniques such 
as direct numerical simulation (DNS) and large eddy 
simulation (LES). DNS is currently limited by com- 
puter hardware to very simple flows at low Reynolds 
numbers, and LES, which directly solves for the large 
turbulent scales and limits empirical modeling to the 
smallest scales, is becoming practical for more com- 
plex flows at higher Reynolds numbers. Birch 4 and 
Bradshaw 5 suggest that LES techniques offer the best 
prospects for improving the capability to calculate tur- 
bulent flows, particularly for flow regions not including 
wall boundary layers. 

As a result, an LES-based technique is an attractive 
option for calculating the mixing dominated regions 
of nozzle flows. However, applying such an LES tech- 
nique simultaneously to the wall bounded regions that 
enter the mixing region (which are an important part 
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of multi-stream nozzles that should he calculated ac- 
curately) will not be practical in the near future. This 
is because computational resources far greater than 
those available today would be required to capture the 
wide range of turbulent time and length scales that, are 
important in such a problem. These turbulent scales 
range from very small eddies in the wall boundary lay- 
ers to very large eddies in the developing mixing layer. 

Nearly all LES simulations of jet and mixing layer 
flows performed to date have placed the inflow of 
the computational domain downstream of any wall 
bounded regions and have either ignored the upstream 
boundary layer effects or used some approximation to 
initialize the turbulent mixing layer. Several authors, 
such as Ragab 6 ' and Hedges 8 have imposed hyper- 
bolic tangent mean velocity profiles at the plane which 
represents the end of t he wall boundary layer regions 
and the beginning of the mixing region. Ragab then 
used the results of a linear stability analysis to generate 
perturbations about the mean velocity profile located 
at the mixing plane. Hedges added small amplitude 
perturbations to the vertical velocity component in 
simulations of heated jet flows. 

The difficulty in using an artificially generated in- 
flow, such as that, assuming a hyperbolic tangent veloc- 
ity profile, is that the characteristics of the upstream 
boundary layers, including velocity, temperature, and 
turbulence profiles, are not accurately represented. 
This is a significant deficiency since the state of the in- 
coming boundary layers have been shown to have sub- 
stantial effects on the development of turbulent mixing 
layers in the experiments conducted by Bradshaw, 9 
Browand and Latigo, 10 and Hussain and Clark. 11 

While RANS- based methods have major deficiencies 
in predicting compressible mixing layers and inher- 
ently are not formulated for calculation of unsteady 
turbulent flows, they have been shown to predict the 
mean flow behavior of wall bounded regions quite well, 
particularly in the absence of adverse pressure gradi- 
ents. As a result, it would be desirable to combine 
a RANS- based technique for the wall boundary layers 
upstream of the mixing region with an LES-based tech- 
nique for the downstream unsteady, turbulent mixing 
region. The development of such a hybrid RANS/LES 
approach is the subject of this work. 

The method developed here is proposed as an al- 
ternative computational technique to performing LES 
calculations everywhere in the computational domain, 
that includes the mean flow characteristics of the in- 
coming boundary layers and is also feasible when con- 
sidering foreseeable computational resources. 

Hybrid RANS-LES Methods 

The realization that LES calculations of flows 
in aerospace and industrial applications at realistic 
Reynolds numbers will not be possible for some time 
has led to interest into the development of hybrid 


techniques. The objective of a hybrid method is to 
retain the essential features of the LES method, but 
to employ a computationally cheaper RANS method in 
regions where it is appropriate. As a result, nearly all 
hybrid methods proposed to date apply a RANS ap- 
proach to attached wall boundary layer regions and an 
LES approach to regions of large scale separation. The 
work detailed in this paper represents the first hybrid 
method development for application to compressible 
mixing layers. 

The most widely publicized hybrid method to date 
is the Detached Eddy Simulation (DES) method of 
Spalart. 12 14 In the DES method, the wall bounded 
regions are calculated using RANS with the Spalart- 
Allmaras 15 one equation turbulence model. Con- 
stant inescu and Squires 10 have applied Spalart's DES 
method to turbulent flow over a sphere, which is an 
appropriate geometry for the method due to the large 
scale separation in the wake of the sphere. 

Speziale 1 * suggested an approach that allows for 
computations varying from RANS in the coarse grid 
limit, through LES, and finally to DNS in the very fine 
grid limit. A Reynolds Stress model is used to close 
the turbulent stresses in the RANS limit, and provides 
the basis for a subgrid model necessary in LES simula- 
tions. Batten et al. 18 also propose a hybrid model that 
employs a Reynolds-Stress model to close the RANS 
and LES equations. Lastly, Arunajatesan et al. 19 have 
applied a hybrid RANS-LES method to cavity flow- 
fields. Their approach employed a two-equation k-kl 
turbulence model to close the RANS equations and 
a one-equation model solving for the filtered subgrid 
kinetic energy to close the LES equations. 

The hybrid RANS-LES method presented here is de- 
veloped for application to configurations such as the 
mixing layer shown in figure 1. This relatively simple 
configuration is representative of more complex nozzle 
geometries in that two wall bounded regions provide 
isolated flows to a single region where compressible 
mixing is the primary flow characteristic. Develop- 
ment of the hybrid method and preliminary assessment 
of the method for a benchmark compressible mixing 
layer configuration are the focus of this paper. 

The hybrid method employs a RANS approach 
to provide the mean flow characteristics of the wall 
boundary layers entering the mixing region. The 
downstream mixing layer is then calculated using LES. 
The method developed here is intended for those noz- 
zle and mixing layer problems in which a geometric 
feature, such as the base region of a nozzle or split- 
ter plate separating the upstream flows, provides the 
dominant unsteady mechanism to drive the develop- 
ment of turbulence in the mixing layer. Although the 
upstream RANS approach does not provide any un- 
steady turbulent information to the mixing layer, the 
mean flow momentum and thermal boundary layer ef- 
fects are calculated and provided to the LES region. 
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Wall Bounded Regions - RANS 


Free Shear Layer Region - LES 


Fig. 1 Schematic of mixing layer demonstrating the hybrid RANS/LES approach 


Method Formulation 

Both the RANS and LES equations are derived 
starting from the general form of the Navier-Stokes 
equations, written in tensor form. The expression 
for conservation of mass, or the continuity equation 
is written in tensor form as: 

^ + dh(pti,-) = 0 (1) 

at oxi 


Conservation of momentum is written: 

d , , d , , dP 

+ ; s; = "as 


If ‘ 2 > 


Conservation of energy is expressed as follows: 


^ + JL( Uj (E t + P)) 


dxj 


9 t \ d 1j IOA 

(3) 


Here, the variable E t represents the total energy (in- 
ternal energy plus kinetic energy) per unit volume: 


Et—pe + -pUiUi (4) 

The equation of state for an ideal gas is used to relate 
the pressure, temperature, and density through: 


P = pRT 


( 5 ) 


For the viscous stresses t x j , it is assumed that the fluid 
is a Newtonian fluid, and as a result, the viscous stress 
is proportional to the rate of strain. This is written: 


where the rate of strain tensor Sij is: 


, _ 1 / dui duj \ 

u 2 dxi ) 


[7) 


Using Stokes's assumption that the thermodynamic 
and mechanical pressures are the same for a fluid un- 
dergoing and expansion or compression: 




(8) 


Equation (6) can then be rewritten as: 

c)u ' 

T ij ~ 2pSij — 3 

The Sutherland model is used to calculate the viscos- 
ity. The heat flux qj is obtained from Fourier's law: 

dT 

* = (10 » 

where k is the thermal conductivity. It is assumed 
that the fluid is thermally perfect, such that the in- 
ternal energy and enthalpy are only functions of the 
temperature, and it is also assumed that the fluid is 
calorically perfect such that the specific heats Cy and 
Cp are constants. As a result, the internal energy and 
enthalpy can be written as: 


e = C V T 
h = C P T 


(11) 


Tij = 2 pSij + A 


duj 

dxj 


Si i 


Assuming that the air is of constant composition and 
does not undergo any chemical reaction, the thermal 
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conductivity i.s only a function of temperature. Using 
the specification of constant specific heats, the follow- 
ing expression is obtained for the thermal conductivity 
as a function of the constant pressure specific heat. 
Prandtl number, and the viscosity: 


_ p£p 

Pi- 


rn 


Mass- Weighted RANS Equations 

In the classical form of Reynolds averaging, the time 
dependent form of the Navier-Stokes equations given 
by (1) through (3) are averaged over a period of time 
that is much larger than the period of turbulent fluc- 
tuations. Each of the dependent variables appearing 
in t hese equations is replaced by the sum of mean and 
fluctuating components. As an example, the velocity 
would be given by: 


Mi = «* -f Uj 


(13) 


where the time averaged velocity u, is given by: 


Mi 



Uidt 


(14) 


For the current work, where fluctuations in density are 
important, a mass (or density) weighting is employed 
in the averaging process, which makes the final form 
of the RANS equat ions much more convenient to work 
with. The dependent variables are again broken into 
mean and fluctuating components: 


Ui = u, + u" (15) 

where the time averaged (using mass weighting) veloc- 
ity lit is given by: 

i r i+T 

Ui = — j puidt (16) 

This mass-weighted Reynolds averaging process is fre- 
quently referred to as Favre averaging, and in general, 
the Favre average of any variable / is defined by: 

f = r=- (17) 

P 

Applying this averaging procedure to equations (1) - 
(3) results in the following expressions for continuity, 
momentum, and energy: 

dp 3 

^ + _ (? „, ) = 0 (18) 


d d dP drij 

dt (pu ' ) + d7 j [pUiUi) + d^ i ~Tx 7 


dr[j 

= ° 


(19) 


I ( Et ) + ^7 ( UjEt + UjP ) - w, {uiTij + UiT > j) 

to + qT P = 0 

(20) 


Spatially Filtered LES Equations 

To derive the LES equations used in this work, the 
time dependent form of the Navier-Stokes equations 
given in equations (1) through (3) is again used as 
the starting point. Instead of time averaging these 
equations, however, an approach similar to the work 
of Ragab and Sheen 6 ' and Erlebacher et, al 20 is used 
that will filter out small scale fluctuations, and only 
retain scales that are large enough to be resolved by 
a particular computational scheme and the computa- 
tional mesh. The filtering operation is defined on any 
variable / by the expression: 

7 (xj)= f G(x-&A)f(&t)<PS ( 21 ) 

Jd 

In equation (21), G is the filter function, D is the flow 
domain, and A is the filter width. The filter width 
A is usually taken to be the grid spacing, and is the 
approach taken here. Note that the overbar used in 
equation (21) indicates a filtered variable. This is in 
contrast to the previous use of an overbar to indicate a 
time averaged quantity in the previous section. As dis- 
cussed by Nelson. 21 the exact form of the filter function 
is not typically known. However, the filter function 
must satisfy: 

/ G{x-tA)d z Z=\ ( 22 ) 

Jd 

In large eddy simulations of compressible flows, it is 
common to use Favre-filtering which is defined as: 

pf 

f= P dr (23) 

P 

where a quantity / is decomposed into resolved and 
unresolved (also referred to as sub-grid scale) compo- 
nents as: 


/ = / + /' (24) 

Equations (23) and (17) are very similar in appear- 
ance, but they refer to very different operations. Both 
operations employ mass (density) weighting, but equa- 
tion (17) defines a time averaged quantity, while equa- 
tion (23) defines a spatially filtered quantity. As was 
the case for the RANS equations, the density, pressure, 
and heat flux terms are not decomposed using mass- 
weighting. Again note that the overline represented a 
time averaging process in the previous section, but it 
will refer to a spatial filtering operation in the current 
section. In addition, Favre filtering differs from Favre 
time averaging in that: 


/#/ 

(25) 

/'# 0 

(26) 
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After applying the filtering procedure to equations (1) 
- (3) the resulting LES expressions for continuity, mo- 
mentum, and energy are: 


dp d ... n 

m + d 7 , (pUi) -° 


(27) 


d dP 

[pui) + — 


dt 


dx; 


dfjj 

dxj 


dxj 


= 0 
(28) 


I (*<) + w s + iiF ) - + 

+!-(?< +o=» 

(29) 


Turbulence Modeling 

Both the RANS and LES sets of equations require 
a turbulence model to close the momentum and en- 
ergy equations. In the RANS approach, all unsteady 
turbulent motion is replaced by a turbulence model. 
The resulting LES equations are very similar in ap- 
pearance to the RANS equations, and also require a 
model to close the momentum and energy equations. 
The difference for the LES equations, however, is that 
the terms replaced by a model are only the turbu- 
lent terms that are too small to be resolved using the 
filtered LES equations. As a result, the large scale tur- 
bulent motion is directly calculated, and the effects of 
the smallest scale turbulence are accounted for using 
a subgrid turbulence model. 

The turbulence model employed here to close the 
RANS equations is the Cebeci-Smith algebraic turbu- 
lence model. 22 23 Since the RANS equations are only 
used in this hybrid method to calculate wall bound- 
ary layer regions with no adverse pressure gradients, 
the selection of a relatively simple algebraic model 
such as the Cebeci-Smith formulation is appropriate. 
The wall function technique of Ota and Goldberg 24 
is used in conjunction with the Cebeci-Smith model 
to enable use of a computational grid with the first 
point off solid boundaries placed in the logarithmic 
layer. This wall function approach is based upon the 
compressible law of the wall formulation of White and 
Christoph. 25,26 The filtered LES equations are closed 
using the Smagorinsky subgrid model. 2 ' 

Implementation of the wall function technique is 
critical to the development of this hybrid approach in 
order to enable use of a single computational grid ex- 
tending continuously from the RANS regions to the 
LES regions. If a wall function approach were not 
used, grids for the RANS regions would have to be 
packed very tightly to the wall and use significant grid 
stretching, while a separate grid which minimizes grid 
stretching would need to be constructed for the LES 


region. Use of such non-continuous grids for the RANS 
and LES regions would require a scheme that would 
likely introduce interpolation errors into the combined 
hybrid method. 

Use of the Cebeci-Smith model to close the RANS 
equations and the Smagorinsky subgrid model to close 
the LES equations is desirable in terms of code im- 
plementation. While the function of the Cebeci-Smith 
model to replace all of the turbulent stresses with a 
model is quite different from that of the Smagorinsky 
subgrid model, which only replaces the small subgrid 
turbulent stresses, both are eddy viscosity models and 
are derived at least in part from mixing-length theory. 
The similar formulation of these two models enables 
the RANS equations and LES equations to be solved 
with a single solution scheme and computational grid, 
as mentioned previously. For a compressible nozzle or 
mixing layer flow, such as that depicted in figure 1, 
the change from RANS regions to LES region occurs 
at the vertical plane passing through the trailing edge 
of the splitter separating the wall bounded flows. 

RANS Turbulence Model 

The unclosed terms from the RANS momentum and 
energy equations are the Reynolds stress and turbulent 
heat flux. The Boussinesq approximation is used to 
relate the turbulent Reynolds stress to the mean rate 
of strain tensor through a turbulent (or eddy) viscosity. 
This is directly analogous to equation (9) which relates 
the viscous stress to the mean rate of strain through 
the molecular (or laminar) viscosity. The turbulent 
analogy to equation (9) is: 


= ti T \2Sii 



(30) 


Similarly, the turbulent heat flux is related to the tem- 
perature gradient through a turbulent conductivity, 


k T : 


</J = pu'jh" 

dT ( 31 ) 

dx j 

The turbulent Prandtl number, Pr T is used to relate 
the turbulent viscosity to the turbulent conductivity: 

Pr r = aL £p (32 ) 

The turbulent Prandtl number is taken to be a con- 
stant here and equal to 0.9. Using equation (32) and 
assuming a constant turbulent Prandtl number enables 
the turbulent heat flux to be expressed as a function 
of the turbulent viscosity that, is used to calculate the 
Reynolds stress. The turbulent heat flux becomes: 


C P p T dT 
Pr T dxj 


(33) 
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The Cebeci-Smith model, which was chosen here to 
close the RANS equations in the wall boundary layer 
regions, treats the wall boundary layer as having in- 
ner and outer regions where the turbulent viscosity is 
defined as: 




T 



y < Vm 

y > ym 


(34) 


In equation (34), y m is defined as the smallest value 
of y (the distance away from a wall) at which fP nner = 
Pouter- The expressions for the inner and outer layer 
turbulent viscosities are as follows: 

Inner Layer: 



with the mixing length t miX is given by: 

fm.x = «;v (l - f -!/+/4+ ) (36) 

Note that the specific form of equation (35) is for a 
two-dimensional boundary layer. 

Outer Layer: 


tauter = 0 tpS v u e Fkieb (37) 

In equation (37), the quantity S v is the velocity thick- 
ness, u e is the boundary layer edge velocity, and F^ieb 
is the Klebanoff intermittency function. The velocity 
thickness is defined as: 



This velocity thickness is identical to the displacement 
thickness for incompressible flows. Klebanoff 28,29 pre- 
sented an expression for the intermittency of turbu- 
lence near the edge of a boundary, which has a func- 
tional form involving the complimentary error func- 
tion. This original intermittency function is usually 
approximated by the following formula, as indicated 
by Cebeci : 22 


f^kleb — 


1 + 5.5 



(39) 


The closure coefficients appearing in equations (36) 
and (37) are 


k = 0.40 a = 0.0168 A+ = 26 

Wall Function Implementation: The Cebeci-Smith 

model is usually integrated down to the wall, using a 
computational grid with the first point off of the wall 
placed well within the laminar sublayer, correspond- 
ing to y + < 5. For the hybrid method developed in 
this work, the objective is to place the first point off 


of the wall in the logarithmic layer to enable the use 
of computational grids that are not packed as tightly 
to the wall. Removing the tight spacing requirement 
will enable a continuous grid into the LES region. In 
addition, because the allowable time step of the com- 
putations is directly proportional to the size of the 
smallest grid cell, a less tightly packed grid enables 
a larger time step for the solution scheme. The wall 
function technique of Ota and Goldberg 24 is one of the 
more simple and effective methods currently in use, 
and it is the technique used in this work. 

Wall functions have been implemented most fre- 
quently in conjunction with two-equation k-e models. 
The benefits of implementing a wall function for use 
with a k-f model are the same as that for the Cebeci- 
Smith model used in this work including reducing grid 
requirements, and increasing the permissible time step 
of the computations. 

The use of a wall-function approach is strictly only- 
valid in flow regions absent of adverse pressure gra- 
dients and separations, due to the assumption that 
the law of the wall holds. However, Avva et al . 30 
have shown results for separated flows in which wall 
function methods perform no worse than methods in- 
tegrating to the wall. The intention of the wall func- 
tion implementation in this work is to only apply the 
method to attached wall boundary layers where the 
law of the wall is valid. 

The Ota-Goldberg wall function employs the White- 
Christoph 25 26 compressible law of the wall: 


sin (— (40) 
yfy V K yt J 

where the compressibility parameter 7 is given by: 


7 = 



(41) 


In equation (40) ut is the value of u + at the first point 
off of the wall, yt is the value of y+ at the first point 
off of the wall, and yt — 0.1287. In equation (41), the 
parameter r is the recovery factor, which is typically 

taken to be Pr3 for turbulent boundary layers, and 
T w is the w r all temperature. An iteration procedure 
is used with equations (40) and (41) to solve for ut , 
from which the shear velocity u T can be obtained: 


u 


+ 

2 


U r t 

U T 


(42) 


Finally, the shear velocity is used to compute the wall 
shear stress through: 


= pu 2 r (43) 

The wall shear stress calculated in equation (43) is 
then used in the solution scheme for the momentum 
and energy equations in the RANS regions. 
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LES Subgrid Scale Model 

The terms that must be closed for LES equations 
are the subgrid-scale stress and the subgrid scale heat 
flux. The earliest subgrid scale model for LES com- 
putations was developed by Smagorinsky. 2 ' Despite 
significant efforts to develop more sophisticated sub- 
grid scale models, the Smagorinsky formulation is still 
widely used, and is itself the foundation upon which 
some of the more sophisticated models are derived. 
The form of the model is very similar to the Cebeci- 
Smith model used for the RANS equations, in that a 
gradient-diffusion mixing-length approach is used. 

The Smagorinsky expression for the subgrid scale 
stress is: 

T’f = p ( U{Uj - UiUj) 

= 2(Cs&fp\/n (Sij - ^SkkSij') - IC'iA-pnSij 

(44) 

The parameter 7 r is defined: 

7T = SijSij (45) 


where the subgrid scale turbulent viscosity is given by: 

/»**' = p(C s A fV* (49) 

Note the similar form of equation (49) to the ex- 
pression for the Cebeci-Smith inner region turbulent 
viscosity in equation (35). While the mixing length 
defined for the Cebeci-Smith model is used to char- 
acterize all of the turbulent motion, the length scale 
defined here for the Smagorinsky model only charac- 
terizes the subgrid-scale motion. Finally, the subgrid 
scale heat flux is modeled analogously to that done for 
the turbulent heat flux of the RANS equations: 



where k sgs is related to /j ags through the turbulent 
Prandtl number. As in the RANS regions, the tur- 
bulent Prandtl number is assumed to be constant in 
the LES regions and equal to 0.9. The subgrid scale 
heat flux becomes 


The parameter A is the filter width and as a result, it 
is also used as the length scale that is characteristic of 
the subgrid turbulence. For use with a computational 
method, A is usually taken to be the grid spacing. In 
three dimensions, with a computational grid having 
unequal spacing in the three directions, this subgrid 
length scale is usually taken to be: 

A = (AxAyAz)^ (46) 

For computational grids with substantially different 
spacing in the three directions, an alternative form (see 
Ragab') is 

\ [(A*) 2 + (At/) 2 + (Azf 

[ 3 

This second expression for the length scale is em- 
ployed in this work. The constants Cs and Cj have 
been found to be dependent on the flow under inves- 
tigation. Rogallo and Moin 31 suggest a range for Cs 
in the range 0.10 < Cs < 0.24. The upper limit on 
Cs given by Rogallo and Moin is investigated for the 
mixing layer in this work. The constant Cj is usually 
equal to 0.01, but several authors, including Ragab' 
and Choi et al. 32 mention that the contribution of the 
term involving Cj may not be important and may be 
neglected. This approach is taken in this work, and as 
a result, the original expression for the Smagorinsky 
subgrid scale stress in equation (44) may be rewritten 
as follows: 

T-J- = 2p' g ’Sij - (48) 


2 

(47) 


Cpji^_ df 
Pr T dxj 


(51) 


Solution Procedure 

The similar form of the RANS and LES equations 
derived here enable both sets to be solved with a single 
computational method. In this work, the Gottlieb- 
Turkel scheme 33 was used. In addition, the equations 
are transformed to generalized coordinates and solved 
using a set of metric terms that are consistent wdth the 
order of the Gottlieb-Turkel scheme. This was done 
to enable the use of the hybrid method on stretched, 
non-Cartesian grids. The Gottlieb-Turkel scheme is 
cited by several authors, such as Hudson and Long, 34 
as providing second order accuracy in time and fourth 
order accuracy in space. Bayliss et al. 35 indicated 
that the scheme has fourth order accuracy only if At 
is of the order (A#) 2 , and Nelson 21 showed that the 
spatial accuracy of the scheme is strictly only third 
order for CFL numbers approaching 1. In the case of 
simulations with grids employing significant stretch- 
ing, however, the spatial accuracy is effectively fourth 
order for most of the computational domain. 


Mixing Layer Simulations 

Experimental Configuration 

The hybrid RANS/LES method has been applied to 
one of the benchmark compressible mixing layer ex- 
periments of Goebel and Dutton 36-38 in which two 
isolated supersonic streams, separated by a splitter 
plate, provide the flows to a constant area mixing 
section. In particular, case 2 of their experiments is 
examined here. A simplified schematic of their ex- 
perimental configuration provided in figure 2 shows 
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top flow bottom flow 


Mach No. 

1.91 

1.36 

c (»'/*■) 

700 

399 

T,(K) 

578 

295 

T(K) 

334 

215 

a(m/.s) 

366 

293 

P(kP a) 

49 

49 

p(kg/m 3 ) 

0.51 

0.79 

S(mm) 

2.9 

2.5 

S* ( mm) 

0.90 

0.44 

9 (mm) 

0.29 

0.21 


Table 1 Flow conditions for case 2 of the Goebel- 
Dutton experiments 

that two isolated streams, in which boundary layers 
develop over a. splitter plate surface, are brought to- 
gether into a constant area mixing section. In all of 
their experiments, the higher speed primary stream oc- 
curred over the top surface of the splitter plate. The 
top stream enters the mixing section axially, while the 
bottom stream enters the mixing section at an angle 
of 2.5 degrees. The splitter plate thickness has a base 
height of 0.5 mm at the trailing edge. Upstream of the 
straight sect ions for the two isolated flows shown in fig- 
ure 2. contoured nozzle blocks provided the supersonic 
flows with nearly uniform exit flow conditions. 

The mixing section height was 48 mm , and the over- 
all length of the mixing section available for flowfield 
measurements was 500 mm. The width of the mixing 
section was 96 mm, and as a result, the mean flow de- 
velopment could be considered two-dimensional. The 
divergence angle of the lower and upper walls of the 
mixing section were adjusted in each experiment with 
two incoming supersonic flows, to account for bound- 
ary layer growth along these two surfaces and to effec- 
tively remove any streamwise pressure gradient. 

Single components LDV measurements w T ere used to 
calculate the boundary layer, displacement, and mo- 
mentum thicknesses of the two streams as they entered 
the mixing section. These quantities are provided 
along with the other operating conditions of case 2 in 
table 1. The documentation of the incoming bound- 
ary layer characteristics makes the Goebel-Dutton ex- 
periments one of the more thoroughly documented 
benchmark data sets available for compressible mixing 
layers. In the mixing region, a two-coinponent LDV 
system was used to measure the axial and transverse 
velocities. In addition, a Schlieren system with a 20 ns 
pulse duration was used to obtain nearly instantaneous 
snapshots of the mixing layer. 

RANS Validation 

To begin the investigation of the hybrid method, 
initial RANS calculations for the two streams feeding 
the mixing layer were constructed. For the Mach 1.91 
stream, wall- integration and wall-function grids were 
generated. Both grids extended 300 mm axially and 


150 mm vertically. In addition, both grids used 141 
points in the axial direction and 141 points in the ver- 
tical direction. The wall-integration grid was packed 
to the wall such that the first grid spacing was 0.006 
mm, corresponding to an average y + of 2.5. The wall 
function grid had the first point placed at 0.05 mm, 
corresponding to an average y~*~ of approximately 20. 
This wall spacing was chosen for use with the wall- 
function grid so t hat the initial grid spacing at the wall 
was exactly l/10th of the splitter plate thickness in the 
experiment. This grid spacing could then be continued 
into the mixing region, with 10 points spaced equally 
in the vertical direction at the base of the splitter. 
In the axial direction, both the wall-integration and 
wall-function grids were packed to leading and trailing 
edges with spacings at the two ends set to 0.10 mm. 
corresponding to 1 /5th of the splitter base thickness. 

For the Mach 1.36 case, grid construction was very 
similar to that used for the Mach 1.91 stream, with the 
exception that the axial dimension was shorter for the 
Mach 1.36 calculations. This was done because the 
Mach 1.36 boundary layer thickness in the Dutton- 
Goebel experiment measured at the beginning of the 
mixing section was smaller than the Mach 1.91 bound- 
ary layer. The grids extended 200 mm in the axial 
direction and 150 mm vertically. Both grids used 141 
points in the axial and vertical directions. The wall- 
integration grid was packed such that the first spacing 
was 0.006 mm, corresponding to an average of ap- 
proximately 3.0. The wall function grid placed the first 
point at 0.05 mm, corresponding to an average y+ of 
approximately 25. 

For each of the cases discussed in this section, the 
wall boundary condition for the calculations was set to 
be an adiabatic no-slip surface and the extreme ver- 
tical boundary was set as a slip surface. The inflow 
boundary was fixed to uniform supersonic flow and 
the outflow boundary extrapolated all quantities. 

Figure 3(a) compares velocity profiles obtained with 
the two approaches to analytical predictions from the 
compressible boundary layer method of Tucker 39 for 
the Mach 1.36 boundary layer. Figure 3(b) provides 
the same comparison for the Mach 1.91 boundary 
layer. The solutions provided with and without the 
wall-function method provide close agreement with the 
Tucker theory. Figure 4 provides a comparison of the 
momentum thickness growth along the two flat plates 
and also indicates close agreement between the cal- 
culations and the method of Tucker. These results 
verify that the RANS part of the current method 
accurately predicts the mean flow behavior of the su- 
personic boundary layers entering the mixing section. 

Two-Dimensional Mixing Layer Calculations 

While true LES simulations require computations in 
three spatial directions, initial two-dimensional hybrid 
calculations of the entire mixing layer configuration 
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Mixing section 


U, =700 m/s 
T m = 578 K 
M 1 = 1.91 


-► 


U 2 = 399 m/s 
T t . 2 = 295 K 
M 2 = 1.36 



Splitter Plate 


A V 


48 mm 


A V 


Splitter Plate Parameters: 

T.E. thickness = 0.50 mm 
Bottom surface angle = 2.5 ° 


not to scale 


Fig. 2 Schematic of Goebel-Dutton mixing layer experiment and Case 2 operating conditions 


were used to investigate effects of axial grid resolution 
and the splitter plate treatment. No subgrid model 
was used in these preliminary calculations. 

The development of the computational model be- 
gan by using the results of the Mach 1.91 and Mach 
1.36 boundary layer simulations. Specifically, the wall- 
function solutions were used, because the computa- 
tional grids employed with the wall-function approach 
enabled a continuous grid into the mixing region for 
use with the hybrid RANS-LES solver. The objective 
was to construct two RANS regions that would pro- 
vide boundary layer quantities S, 6*, and 0 that nearly 
matched those measured in the experiment. Because 
it would be virtually impossible to match all of three 
quantities exactly, the momentum thickness 0 was cho- 
sen as the key parameter to match the computations 
with the experiment. The momentum thickness repre- 
sents the mean momentum deficit entering the mixing 
section and is fundamental to the mixing layer behav- 
ior. 

Examining the Mach 1.91 boundary layer first, ta- 
ble 1 indicates that the momentum thickness for this 
stream was measured as 0.29 mm at the trailing edge 
of the splitter plate. For the wall function calcula- 
tions discussed in the previous section, the momentum 
thickness reached 0.29 mm at a Reynolds number of 
3,680,000 (see figure 4(b)), corresponding to an axial 
position of 198 mm from the leading edge of the plate. 

The next step was to construct a new computational 
grid. The axial domain was shortened to 198 mm while 
retaining the same number (141) of axial grid points. 
The grid stretching was modified to accommodate the 
shorter domain while maintaining the initial and ter- 


minal grid spacings. The vertical domain was reduced 
to 23.75 mm to exactly match the height of the Mach 
1.91 stream in the experiment. The vertical domain 
of the original grid reached 23.75 mm at the 94th grid 
point, so the first 94 points from the original grid were 
used in the modified grid. Calculations obtained with 
this 141 x 94 point, 198 mm by 23.75 mm grid pro- 
vided boundary layer quantities identical to that of 
the original 141 x 141 grid, further validating the grid 
independent characteristics of the RANS method. 

A similar procedure was used for the Mach 1.36 
boundary layer. Table 1 indicates that the momentum 
thickness for this stream was 0.21 mm at the split- 
ter trailing edge. Examination of the wall-function 
solution obtained for the Mach 1.36 boundary layer re- 
vealed that the momentum thickness became 0.21 mm 
at an axial position of approximately 120 mm, corre- 
sponding to a plate Reynolds number of 2,720,000. A 
modified grid was constructed using 141 axial and 94 
vertical points for this Mach 1.36 case, corresponding 
to a physical domain of 120 mm by 23.75 mm. Calcu- 
lations with this grid provided a solution identical to 
that obtained with the original 141 x 141 grid. 

The last step before constructing the entire RANS- 
LES computational grid was to extract the grid points 
from the modified wall-function grids just discussed, 
extending from 81-141 in the axial domain. The solu- 
tions were then held constant at the 81st axial station 
as the inflow of the final hybrid RANS-LES grid. This 
was done to reduce the grid requirements for the subse- 
quent RANS-LES computations. Boundary layer cal- 
culations subsequently obtained with these two short- 
ened axial grids using 61 axial points and 94 vertical 
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Fig. 3 Velocity profiles for boundary layer calcu- 
lations 


Fig. 4 Momentum thickness growth along super- 
sonic flat plates 


points again returned the same boundary layer quanti- 
ties and velocity profiles as the other solutions. These 
were the RANS grids that were used to join with the 
LES mixing region computational domain. 

Three axial grid spacings (with 200, 400, and 800 
points in the mixing section respectively) were exam- 
ined in the initial two-dimensional hybrid calculations. 
In all cases, however, the vertical spacing from the 
two RANS regions was continued throughout the en- 
tire LES region. With the tightest vertical spacing of 
the wall-function boundary layer solutions set to 0.05 
mm, 10 grid spacings were used vertically across the 
splitter base. As a result, all of the hybrid grids used 
197 vertical points in the mixing region. 

Fixed inflow boundary conditions were used at the 
RANS inflows. The fixed inflow for the Mach 1.91 
upper stream was placed at an axial position 67 mm 
upstream of the splitter plate trailing edge and the 
Mach 1.36 lower stream inflow was placed 42 mm up- 


stream of splitter tip. At the outflow of the mixing 
section, corresponding to an axial distance of 300 mm 
from the trailing edge of the splitter plate, an extrap- 
olation boundary condition was used, which is appro- 
priate for the mixing supersonic flow exiting the axial 
domain. The top and bottom walls of the mixing sec- 
tion were approximated as slip walls, and no attempts 
were made to simulate boundary layers developing on 
these two surfaces. The Goebel-Dutton experiment 
configuration was specifically designed with adjustable 
divergence angles for these walls to account for bound- 
ary layer growth, and to thereby provide data that may 
be directly compared to calculations which do not in- 
clude the mixing section boundary layers. 

The first two-dimensional simulation investigating 
grid density effects was for the coarsest grid using 200 
axial points. A vortex shedding pattern was observed 
to originate from the trailing edge of the splitter plate 
due to the separation of the two flows leaving the wall 
bounded RANS regions and entering the LES mixing 
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section. The vortex shedding quickly dissipates and 
the flow appeared to be mostly laminar throughout 
the rest of the mixing section. 

The second two-dimensional simulation was for the 
computational grid using 400 axial points. The char- 
acteristics of t his flow were substantially different from 
those of the 200 grid point case. Instantaneous density 
contours are provided in figure 5. A stronger vortex 
shedding is evident for this case, although the vortex 
strength gradually dissipates back to a laminar state 
approximately 40 mm downstream of the splitter trail- 
ing edge. An instability forms at an axial position of 
80 mm which in turn slightly dissipates before resum- 
ing growth at a position 180 mm downstream of the 
splitter plate trailing edge. A turbulent pattern then 
grows from this location to the exit at x =300 mm. 

The 400 point axial grid consisted of smaller axial 
grid steps than the 200 axial point grid and had a 
significantly reduced stretching factor relative to the 
coarse grid. As a result, the truncation error term in 
the Gottlieb-Turkel scheme that effectively smooths 
discontinuities is substantially reduced for the 400 ax- 
ial point grid, and the capability to resolve unsteady 
flow behavior was improved. 

The third computational grid investigated had 800 
axial points. A substantially different flow behavior 
was also observed for this case compared to the so- 
lutions obtained with 200 and 400 axial points. The 
density contours in figure 6 again indicate a vortex 
shedding pattern that originates from the trailing edge 
of the splitter plate, but unlike the 400 point case, the 
solution does not return to a laminar state before tran- 
sitioning over to a turbulent-like pattern at x = 80 
mm. The very tight axial spacing for this case is suffi- 
cient to minimize the truncation error damping effects 
on the unsteady flow development. 

An additional two-dimensional grid with 800 axial 
points was constructed for a modified splitter geometry 
in which the splitter trailing edge is reduced to a sharp 
tip, and as a result, the flow separation and vortex 
shedding is removed. Figure 7 provide instantaneous 
density contours for the case with a sharp trailing edge. 
As expected, the vortex shedding evident in the base- 
line case was removed in the current case with the 
sharp tip. The lack of a separation region results in an 
initially laminar mixing layer through the beginning 
of the mixing section. Interestingly, the laminar flow 
begins to transition to a more turbulent structure at 
nearly the same position observed for the baseline case, 
at approximately x = 80 mm. The structures remain 
relatively small until x — 150 mm where large scale 
turbulence forms. These structures are more similar 
to the Brown-Roshko organized structures than were 
those of the baseline case with 800 axial points. 


Three-Dimensional Mixing Layer Calculations 

In the previous section, two-dimensional calcula- 
tions were used to construct the initial computational 
model of the mixing layer and to examine preliminary 
effects of grid resolution and splitter plate treatment. 
To correctly investigate the capability of the hybrid 
method, however. LES calculations obtained in three 
spatial directions with the use of a subgrid scale model 
were performed. 

The grid topologies and boundary conditions used 
for the three-dimensional simulations were very simi- 
lar to those used for the two-dimensional simulations. 
The two-dimensional computational grids with 200, 
400, and 800 axial point grids were used to con- 
struct the three-dimensional grids used here. To add 
the third computational direction in each case, the 
two-dimensional planar grid was copied to provide 11 
points in the third (or z) direction. The grid spac- 
ing in the z direction was uniform and set equal to 
the axial spacing at the splitter trailing edge, or Az = 
Ax i = 0.10 mm. Because of the very small number of 
grid points used in the z direction and the small phys- 
ical space that is represented, only very small wave 
components in this direction could be simulated, and 
a periodic boundary condition was used in this direc- 
tion. 

All of the other boundary conditions and the so- 
lution procedure are identical to that used for the 
two-dimensional mixing layer calculations, with the 
one exception being that the Smagorinsky subgrid 
scale model was used in these three-dimensional sim- 
ulations. The switch from the RANS regions to the 
LES region at the mixing plane (corresponding to a 
vertical plane drawn through the trailing edge of the 
splitter plate) was accomplished by changing the eddy 
viscosity used in the flow solver from the Cebeci-Smith 
turbulence model to the Smagorinsky subgrid scale 
model. As a result, the effect of the eddy viscosity 
changes from that of replacing all of the turbulent 
stresses in the RANS regions to that of only replac- 
ing the subgrid stresses in the LES regions. 

The first three-dimensional simulation was obtained 
using the computational grid with 200 axial points. 
The behavior of the flow just downstream of the split- 
ter trailing edge is very similar to the corresponding 
two-dimensional calculation. The lack of adequate ax- 
ial grid resolution results in a rapid dissipation of the 
initial vortex pattern, and is even more rapid in the 
three dimensional case due to the dissipative nature of 
the Smagorinsky subgrid scale model. 

The three-dimensional grid using 400 axial points 
was utilized next. The density contours shown in fig- 
ure 8 indicate a fundamentally different flow structure 
than observed with any of the two-dimensional calcu- 
lations. In particular, the vortex shedding is observed 
to transition into a turbulent pattern much closer to 
the trailing edge. While the Schlieren photographs 


N AS ATM- 200 1 -2 1 0762 


11 




a) Beginning of mixing section 



b) Entire mixing section 

Fig. 5 Instantaneous density contours for the 400 axial grid point case (2D) 



a) Beginning of mixing section 



b) Entire mixing section 

Fig. 6 Instantaneous density contours for the baseline 800 axial grid point case (2D) 
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a) Beginning of mixing section 



b) Entire mixing section 


Fig. 7 Instantaneous density contours for the 800 axial grid point case with a sharp trailing edge for the 
splitter (2D) 



a) Beginning of mixing section 



b) Entire mixing section 

Fig. 8 Instantaneous density contours for the 3D case using 400 axial points 
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Fig. 9 Velocity profiles in the mixing section 
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from the Goebel- Dutton experiments did not probe the 
flow details near the splitter plate trailing edge, such 
details were examined by Clemens and Mungal 40 for 
a similar experiment configuration and flow operating 
conditions. A Schlieren photograph from the Clemens- 
Mungal experiments indicates an initial flow structure 
similar to that observed in the three-dimensional cal- 
culation with initial vortex shedding from the trailing 
edge of a splitter plate followed by a rapid transition 
into turbulence. The turbulent structure in the ex- 
periment was observed to be of primarily small scales, 
while the LES calculations, by definition, only cap- 
ture the large scale structures. While the length of 
the organized vortex structure in the calculations does 
not exactly match that of the Clemens-Mungal ex- 
periment, the three-dimensional calculations predict 
the transition location much more accurately than the 
previous two-dimensional results. In addition, these 
results verify that LES calculations must be run in 
three dimensions to allow instabilities to form in all 
three directions. 

Mean axial profiles obtained from this three- 
dimensional calculation are compared to experimental 
data of Goebel and Dutton in figure 9. The calculation 
indicates greater shear layer spreading than shown by 
the experimental data. The profiles of u r ms in figure 
10(a) and v rms in figure 10(b) generally indicate over- 
predictions from the calculation, which corresponds 
to the wider axial velocity profiles in figure 9. Liou 
et al. 41 and Inoue 42 also reported large overpredic- 
tions in the turbulence intensities for planar mixing 
layers. With the confinement of the current three- 
dimensional calculations to a very small domain in 
the z direction, the inability to calculate large scale 
fluctuations in this direction may be responsible for 
the overpredictions of u rms and v rm$ . In addition, the 
Goebel-Dutton Schlieren photographs indicated a very 
fine turbulent structure contained within the larger 
scale development, while the LES calculations inher- 



Fig. 10 Turbulence intensities in the mixing sec- 
tion 

ently can resolve only the larger turbulent scales. 

The last three-dimensional case was run using 800 
axial points in the mixing section. The finer grid used 
reduced the permissible time step by nearly a factor of 
two relative to the 400 axial point case. This reduced 
time step and doubling the number of grid points in the 
mixing section would require a factor of four increase 
in the computer CPU time requirements to run this 
case to completion, relative to the case with 400 axial 
points. Considering that the 400 axial point case re- 
quired 500 CPU hours on a Cray C90 computer, 2000 
Cray C90 hours would be required to complete the 
case with 800 axial grid points. As a result, this last 
three-dimensional case was run long enough to allow 
the flow to fully develop in the mixing section, but not 
long enough to enable time averaging of the turbulent 
statistics. 

In contrast to the grid refinement studies performed 
for the two-dimensional cases, the large scale turbulent 
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development did not change significantly when increas- 
ing the number of axial points in the three-dimensional 
computations from 400 to 800. The instantaneous 
density contours in figures 11 indicate a large scale 
turbulent structure which closely resembles those of 
the 400 axial grid point case. In particular, the break- 
down of the organized vortex structure to turbulence 
is very similar to those previously shown for the 400 
axial grid point case. Within the large scale struc- 
tures, more fine scale turbulence is evident in the 800 
axial point case. This behavior corresponds directly 
with the philosophy of LES in that as the computa- 
tional grid is refined, smaller structures are able to 
be resolved and the role of the subgrid scale model is 
reduced. In the idealized limit of a grid which is suf- 
ficiently fine to resolve all turbulent length scales, a 
direct numerical simulation is obtained. 

Conclusions 

The work described here represents the initial efforts 
to develop and evaluate a hybrid RANS-LES method 
for compressible mixing layer simulations. Although 
the RANS approach does not provide any unsteady 
turbulent information to the LES region, the mean 
flow boundary layer characteristics are provided. The 
hybrid method was developed for the analysis of nozzle 
and mixing layer configurations in which a dominant 
structural feature, such as the base region of a nozzle 
or splitter plate separating the upstream flows, will 
provide the dominant unsteady mechanism to drive 
the development of turbulence in the mixing layer. 

The Cebeci-Smith turbulence model, despite its rel- 
atively simple form, was demonstrated to provide ac- 
curate calculations of boundary layer flows in the 
RANS regions that are free of adverse pressure gra- 
dients or separations. Further, the use of the Cebeci- 
Smith model in conjunction with the Ota-Goldberg 
wall function enabled calculations of supersonic wall 
boundary layers to nearly the same accuracy as that of 
the standard approach of integrating the Cebeci-Smith 
model through the viscous sublayer, while enabling a 
significantly larger vertical grid spacing near the wall. 
The wall function approach enabled a continuous com- 
putational grid to be used from the RANS to the LES 
regions, and the method thereby avoided the use of 
discontinuous grid zones that would have required an 
interpolation scheme at the RANS-LES interface. In 
addition, the origins of the Cebeci-Smith RANS tur- 
bulence model and the Smagorinsky LES subgrid scale 
model are both in mixing length theory, and this simi- 
lar form of the models assisted in code implementation. 
As a result, the use of a more sophisticated turbulence 
model to close the RANS equations was found to be 
unnecessary, provided the RANS regions are restricted 
to attached, zero pressure gradient wall boundary layer 
regions. 

While true LES calculations require computations in 


three spatial directions, two-dimensional simulations 
of a benchmark mixing layer experiment were consid- 
ered first to address effects of axial grid resolution and 
splitter plate treatment. The parametric study of ax- 
ial grid resolution indicated more realistic turbulent 
development with increasing axial grid density. For all 
of the cases examined, a vortex shedding was found 
to originate from the base region of a splitter plate 
separating the upstream wall bounded regions. For 
the finest grid examined, the unsteady vortex pattern 
eventually transitioned to a turbulent structure. The 
location of this transition, however, was much fur- 
ther downstream than observed in the experiments. 
Calculations obtained for the case in which the finite 
thickness splitter base was changed to a sharp tip in- 
dicated that the vortex shedding was removed. 

Three-dimensional calculations were obtained next 
for grids constructed by copying the two-dimensional 
planar grids to locations in the third computational 
direction. Only a small domain w'as modeled in 
this third direction, and periodic boundary conditions 
were employed along the extreme boundaries. For 
the coarse three-dimensional grid, again no turbulent 
flow development was observed. For the intermedi- 
ate grid, the vortex shedding found previously in the 
two-dimensional simulations was also observed in the 
three-dimensional calculations. However, the orga- 
nized vortical structure rapidly disintegrated into a 
significantly more realistic turbulent flow structure. 
This rapid transition to turbulent flow was nearly iden- 
tical to that found in experimental investigations of a 
similar mixing layer configuration. Although the ex- 
tent of the third dimension in these calculations was 
very small, an unsteady mechanism by which distur- 
bances could develop in all three directions and result 
in a rapid transition to turbulent flow was enabled by 
the three dimensional calculations. In contrast, a two- 
dimensional approach, by definition, does not allow for 
such three-dimensional disturbances to develop. The 
results of these calculations verified that LES simula- 
tions must be performed in three dimensions. A final 
three-dimensional calculation was investigated using a 
computational grid constructed from the most densely 
packed two-dimensional grid. Because a prohibitively 
long run time would be required to complete this 
solution for turbulent statistics purposes, the calcu- 
lations were run only long enough to allow the flow 
to fully develop. The initial transition to turbulence 
and large scale turbulent structures evident for this 
case were very similar to those for the intermediate 
three-dimensional case. More resolution of the finer 
turbulent scales contained with the larger structures 
was observed for the fine grid case, which is in line 
with the philosophy of LES to resolve finer scales as 
the grid density is increased. 

Comparisons of time-averaged axial velocities and 
turbulence intensities from the calculations to experi- 
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a) Beginning of mixing section 



b) Entire mixing section 

Fig. 11 Instantaneous density contours for the 3D case using 800 axial points 


mental data indicated reasonable agreement, with the 
solutions indicating somewhat higher levels of turbu- 
lent mixing. A major source of discrepancy between 
the calculations and experiment is believed to be the 
lack of adequate grid resolution to resolve the small 
turbulent scales contained within the larger turbulent 
structures. Anot her source of discrepancy was the very 
small domain used in the third computational direc- 
tion. Despite these limitations, the three-dimensional 
calculations demonstrated the success of the hybrid 
method to capture the dominant characteristics of the 
mixing layer, and in particular, the rapid transition of 
the organized vortex structure to a turbulent mixing 
layer structure. It is expected that improvements in 
the fidelity of the solution scheme, and more impor- 
tantly, improvements in computing power, will enable 
better predictions of the turbulent statistics. 
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